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A software program and associated methodology to study gust loading on aircraft exists 
for a classification of geometrically simplified flexible configurations. This program consists 
of a simple aircraft response model with two rigid and three flexible symmetric degrees of 
freedom and allows for the calculation of various airplane responses due to a discrete one- 
minus-cosine gust as well as continuous turbulence. Simplifications, assumptions, and 
opportunities for potential improvements pertaining to the existing software program are 
first identified, then a revised version of the original software tool is developed with 
improved methodology to include more complex geometries, additional excitation cases, and 
output data so as to provide a more useful and accurate tool for gust load analysis. Revisions 
are made in the categories of aircraft geometry, computation of aerodynamic forces and 
moments, and implementation of horizontal tail mode shapes. In order to improve the 
original software program to enhance usefulness, a wing control surface and horizontal tail 
control surface is added, an extended application of the discrete one-minus-cosine gust input 
is employed, a supplemental continuous turbulence spectrum is implemented, and a 
capability to animate the total vehicle deformation response to gust inputs in included. These 
revisions and enhancements are implemented and an analysis of the results is used to 
validate the modifications. 


Nomenclature 


c 

= 

generalized damping matrix 

c 

= 

chord, ft 

K 

= 

generalized stiffness matrix 

L 

= 

scale of turbulence, ft 

M 

= 

generalized mass matrix 

Q 

= 

generalized force matrix 

T 

= 

kinetic energy 

U 

= 

elastic energy 

V 

= 

velocity, m/s 

w g 

= 

gust input, m/s 

a 

= 

angle of attack, radians 

s 

= 

generalized coordinate 


I. Introduction 

C ontinuous efforts to advance technologies and design more sophisticated aircraft have been made since Wilbur 
and Orville Wright's first successful flight in 1903. With such advancements, aircraft evolved into larger, more 
flexible, increasingly complex vehicles that could travel farther, at higher altitudes, and at faster speeds. As the 
evolution in aircraft design progressed, the need arose to also improve techniques for analyzing the characteristic 
behaviors of aircraft. One significant area of consideration is given to the response of aircraft to the turbulent 
atmospheric flight environment. Within this area of study, consideration is given to both the response of the aircraft 
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to a discrete gust encounter and to the response of the aircraft to continuous turbulence. As aircraft design evolved, 
the methods utilized to study the behavior of aircraft in turbulent air also progressed. 

With the engineering goal of designing and operating aircraft safely, reliable methods of determining gust loads 
are required. Since the first National Advisory Committee for Aeronautics (NACA) publication relating to gusts 
loads was written in 1915, the scope of gust analysis procedures has progressed from discrete time -domain analysis 
methods to continuous frequency-domain analysis via Power-Spectral Density (PSD) methods. In chronological 
order, the following gust analysis methods were introduced: 1) Sharp-Edged Gust Formula; 2) Ramp Gradient Gust 
Formula; 3) Triangular Gradient Shape Gust; 4) Revised Gust Formula (Pratt's Method); and 5) Power-Spectral 
Density Analysis. 1 Pratt's method assumes a one-minus-cosine discrete gust profile while the PSD approach assumes 
that the turbulence is a Gaussian random process and allows for study over the entire frequency range for which 
response values are needed while covering multiple degrees of freedom. While the FAA requirements include only 
those associated with continuous turbulence, discrete gust load analysis methods are still considered useful in 
providing insight to aircraft behaviors in response to a turbulent atmosphere and may be used as an alternate means 
of compliance. 

In a 1997 publication titled "A MATFAB Program to Study Gust Foading on a Simple Aircraft Model," a 
program to study aircraft loading due to discrete gusts and continuous turbulence was presented. 2 While this 
software program enables the calculation of five output load quantities due to discrete gusts as well as continuous 
turbulence, the five degree-of-freedom aircraft model is limited to simple geometries. The following geometric and 
structural parameter variations are permitted: 1) aircraft wing area and horizontal tail area; 2) aircraft wing sweep 
angle; 3) tail length; 4) aircraft center of gravity position; and 5) aircraft main mass properties. While the aircraft 
wing area is variable, the restriction exists that the wing root chord and wing tip chord must be equal. For reference, 
the structural model is shown in Fig. 1. Note in Fig. 1 that the aircraft configuration is defined by five invariable 
mass points along the swept-back aircraft wing, ten invariable fuselage mass points, and one invariable mass point 
on the unswept, untapered horizontal tail, limiting the complexity of the structural configuration and unrealistically 
limiting the possible planform geometries. 

In order to improve the original program to enhance its usefulness, a capability to include more complex, user- 
defined geometries is included, a wing control surface and horizontal tail control surface are added, an extended 
application of the discrete one-minus-cosine gust input is employed, a supplemental continuous turbulence power 
spectrum is implemented, and a capability to animate the total vehicle physical deformation response to gust inputs 
is provided. 

The paper begins with a brief description of the original software program and its associated methodology. 
Improvements and enhancements are then discussed, followed by results and concluding remarks. 

II. Original Gust Response Model Program 

A. Aircraft Model Configuration and Mode Shapes 

The original two-dimensional half-aircraft response model consists of a beam-model fuselage, a swept, 
rectangular wing, and an unswept, rectangular horizontal tail." The original structural model is defined by five wing 
mass points, one tail mass point, and ten fuselage mass points with a right-handed axis system having its origin 
situated at the intersection of the spanwise extension to the leading edge of the Mean Aerodynamic Chord (MAC) 
and the aircraft centerline. The global x-axis is aligned with the vehicle centerline and is defined as positive forward, 
the global y-axis is defined to be positive moving in the outboard direction along the right wing, and the global z- 
axis is defined as positive downward. The original aircraft structural model, including the position of the center of 
gravity, is shown in Fig. 1. A local stability axis system lies along the wing elastic axis and is denoted as x' and y' in 
Fig. 1. 

Two rigid body modes and three flexible modes are included in the gust response model program. The rigid body 
displacement modes are defined by a translational motion in the z-direction, known as plunge, and a rotation about 
the y-axis, known as pitch. Translational motion in the x-direction is not considered. The plunge mode, mode 1, is 
positive in the direction of the positive z-axis and the pitch mode, mode 2, is defined by a positive nose up rotation 
about the y-axis. 

The flexibility of the airplane is assumed to be sufficiently well captured by two wing modes (1 st bending and 1 st 
torsion) and one rear fuselage mode (1 st bending). Rear fuselage bending and wing bending are positive in the 
direction of the positive z-axis, and wing torsion is positive leading edge up. The elastic deformation modes are 
given as simple mathematical expressions describing the deformation curve rather than being associated with a set 
of eigenmodes. Illustrations of one rigid body mode and one flexible mode are given in Figs. 2 and 3, respectively. 
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Figure 1. Original Configuration Structural Model 
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Figure 3. Wing Bending Deformation Mode 


B. Mathematical Approach, Inputs, and Output Load Responses 

A convenient approach to determining the equations of motion for systems with multiple degrees of freedom is 
available due to J.L. Lagrange. ’ Derivation of the Lagrange equation is not presented in this paper; however, an 
explicit derivation of this equation is presented in Refs. 3 and 4. For many problems in aeroelasticity, the use of 
Lagrange's equation reduces the complexity of the mathematics involved in deriving the equations of motion by 
eliminating the dependent Cartesian coordinates and introducing independent variables known as generalized 
coordinates. 4 The generalized coordinates implicitly contain the constraints on the motion of the system. By 
transformation of the variables from a Cartesian coordinate system to generalized coordinates, the equations of 
motion can be transformed into the form of the Lagrange equation. The Lagrange equation is expressed as 

d dT 
dt d^i 

For simplicity, time-domain solution procedures are carried out in the frequency-domain via Fourier transforms 
rather than in the time-domain with direct integration; therefore, this program is only appropriate for studying the 
responses of linear systems. Processing of the Lagrange equation results in a set of n second order governing 
equations of motion of the form 

Nll + cl + Ki = Q r Z + Q w w g ( 2 ) 

Solution of the equations of motion expressed in generalized coordinates becomes 

I = A~ 1 Q w w g (3) 

The aerodynamic forces and moments are computed using simple strip theory. In the original program, the 
aircraft wing is divided into five equal strips with a constant wing lift curve slope along the wing span and the 
horizontal tail is represented by one strip with a constant tail lift curve slope along the tail span. The unsteady 
aerodynamic force acting on a thin airfoil undergoing unsteady motion was obtained by Herbert Wagner, Hans 
Kiissner, Theodore Theodorsen, William Sears, and others. 4 In this work unsteady effects are accounted for by 
applying Theodorsen’s function to the aircraft response-induced angle of attack, and by applying Sears' function to 
the gust induced angle of attack. A cogent derivation of these functions is given in Ref. 4. 


dT dU „ 

■ 1 — Q: 

dti 


( 1 ) 
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Excitation approaches include a discrete one-minus-cosine shaped gust and the von Karman continuous 
turbulence power spectrum. An illustration of the one-minus-cosine shaped gust input is shown in Fig. 4. The five 
available output load quantities are: (1) center of gravity acceleration (load factor); (2) wing root shear force; (3) 
wing root bending moment; (4) wing root torsion moment; and (5) tail root shear force. Assuming that each output 
load can be expressed as a linear function of the input load and generalized coordinates, the output load can be 
written as 


Substituting Eq. (3) into Eq. (4) yields 
Setting 

the output load equation becomes 


y = C!;+Dw g 

y = (ca~ 7 q w + d \ 8 

H = CA~‘Q w + D 
y = Hw g 


(4) 

(5) 

(6) 
(7) 


where H is the transfer function. Given H , the outputs due to the inputs can now be determined. An example of 
the five output time histories due to a discrete one-minus-cosine gust input is given in Fig. 4. All responses 
correspond to the structural configuration shown in Fig. 1 with default flight conditions. The default aircraft 
dimensions, main mass properties, and flight conditions for this configuration can be found in Ref. 2. The time 
histories in Fig. 4 exhibit three important characteristics. First, for all quantities, the initial response correlated very 
closely with the 0.4 second duration and shape of the one-minus-cosine gust input. Second, for all quantities, 
following the initial response, the presence of a low frequency response corresponding to the vehicle short-period 
mode is observed. Third, visible primarily in wing shear force and bending moment, following the initial response 
and superimposed on the low-frequency response, the presence of a higher frequency response, corresponding to the 
wing bending mode, is observed. 



Time [s] 


Figure 4. Time Responses due to One-Minus-Cosine Discrete Gust Input 
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III. Program Methodology Revisions and Enhancements 

The objective of this work was to modify the original program developed in Ref. 2 to include enhancements that 
improve it's utility. A brief discussion of the modification procedures and implemented enhancements is given in 
this section. 

A. Modification of Aircraft Structural Model 

The expansion of the existing structural parameters to include implementation of a variable number of mass 
points on the aircraft wing, fuselage, and horizontal tail, additional options for mass, inertia, and bending and 
torsional stiffness distributions, implementation of a variable aircraft wing and horizontal tail taper ratio, and 
implementation of a variable horizontal tail sweep angle is discussed in this section. Additionally, the revised 
structural model is developed in a modified axis system. These modifications improve the gust response model 
presented in Ref. 2 by creating additional user-defined parameters and by eliminating certain program limitations to 
allow for the inclusion of more complex geometries. 

Development of the aircraft model to include more complex geometries requires the capability for 
implementation of a swept-back, tapered aircraft wing and a swept-back, tapered horizontal tail. This feature has 
been added to allow for user-defined sweep angles and taper ratios. In order to develop a more precise structural 
model for analysis of gust loads, it is useful to be able to vary the number of mass points along the span of the 
aircraft wing and horizontal tail, and along the length of the fuselage. Recall that computation of the aerodynamic 
forces and moments is accomplished utilizing simple strip theory aerodynamics; therefore, because there is one 
lumped mass per designated strip i , n desired mass points equates to n number of strips chosen by the user for the 
aerodynamic model. When executing the program, the user is prompted to input the desired number of strips for the 
model, thus simultaneously defining n number of mass points along the aircraft wing, horizontal tail, and fuselage. 
The location of mass lump i is consistent for each strip in the three fundamental elements comprising the aircraft 
model. All strips are of equal width along the span of the respective component; consequently, each lumped mass is 
equally spaced along the span of the aircraft wing and horizontal tail and along the length of the fuselage. A revised 
structural aircraft model example comprised of twenty aircraft wing, horizontal tail, and fuselage lumped masses is 
shown in Fig. 5 in a new global axis system. The illustrated configuration mentioned above contains a planform 
identical to Fig. 1, but employs the previously described modifications, thus creating a structural model that is more 
representative of the vehicle. 

In addition to the structural model modifications, a useful feature is to provide options for the manner in which 
masses and inertias are distributed along the aircraft components. The following distribution options are available to 
the user: 1) default linear distribution; 2) user-defined mass and inertia values per strip; and 3) user-defined 
distribution factors. The default distribution is linear and default distribution factors are used to compute the mass 
and inertia values for each respective strip. The second distribution option allows the user to manually define mass 
and inertia values for each strip. The third option allows the user to define desired distribution factors for each strip. 
From the user-defined distribution factors, given the total component mass and inertia values, the respective values 
for mass and inertia are computed in the same manner as for option 1 . 

Similar to the added features for lumped mass and inertia distributions, bending stiffness, El, and torsional 
stiffness, GJ, distribution options are also implemented. A wing element has a constant bending stiffness and a 
constant torsional stiffness, a fuselage element has only a constant bending stiffness and is assumed to be torsionally 
rigid, and a tail element is considered to be rigid. After testing, the default linear distribution for fuselage bending 
stiffness was found to be reasonable; therefore, the fuselage bending stiffness distribution options follow those given 
for the lumped mass and inertia distributions. While a linear stiffness distribution provides an adequate 
approximation for the aircraft wing, this distribution does not yield the most reasonable fit. Therefore, an additional, 
more accurate approximation is employed. This alternate option utilizes a cubic spline interpolation method to 
determine the stiffness distribution along the aircraft wing. A requirement for this option is that the user defines El 
and GJ values for lumped masses at 10%, 30%, 50%, 70%, and 90% semi-span. Interpolation is then performed to 
obtain the El and GJ values at the remaining lumped mass locations. An example wing bending stiffness distribution 
resulting from the cubic spline interpolation method given an aircraft wing comprised of ten lumped mass points is 
given in Fig. 6. This distribution is plotted over the original distribution found in Ref. 2 and illustrates accurate 
agreement between the two implementations. 

For user convenience, a new global axis system was desired. This new axis system is right-handed and has its 
positive x-axis pointing aft and its positive z-axis pointing up. The origin of this new axis system is located at the 
intersection of the wing leading edge with the fuselage centerline. Therefore, for all straight and aft-swept wings 
(and even for some forward-swept wings), the origin of the new axis system is forward of the origin of the old axis 
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system. Expressed in the original global axis system, the longitudinal offset of the origin of the new global axis 
system from the origin of the old global axis system is 

*o = yMAC tmA ( 8 ) 


where y MAC is the y-location of the leading edge of the MAC and A is the sweep angle of the leading edge of the 
wing. The transformation equations are given as 
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B. Modification of Aerodynamic Strip Model 

It is beneficial to permit the user to define an aerodynamic model comprised of a variable number of strips. The 
procedure implemented in the development of the revised aerodynamic strip model generates an aerodynamic model 
comprised of a user-defined equal number of strips along the aircraft wing, horizontal tail, and fuselage. An example 
illustration of a fifty-strip configuration is provided in Fig. 7. As a result of the improvement applied to the 
generation of the aerodynamic model, it is necessary to modify the respective lift and moment equations to reflect 
the user-defined number of strips. Additionally, whereas the average chord is equal to the MAC in the original 
configuration, the average chord varies per strip of the modified tapered configuration. For a wing with n strips 
there are n + 1 boundaries. For the i' h strip, the inboard and outboard boundaries have chord lengths c, and c i+1 , 
respectively. Therefore, the average chord for the i' h strip is given by 


c 



(10) 


and the lift and moment equations are then solved using the average chord per strip. Computation of the horizontal 
tail average chord per strip follows Eq. (10). Likewise, computation of the lift force acting on the horizontal tail is 
executed using the average chord per tail strip. A difference exists when accounting for the decrease in horizontal 
tail angle of attack produced by the downwash angle. In the original configuration, the downwash angle, was 
computed with reference to the second wing strip. For the five-strip original configuration, the y-location of the 
outboard boundary of the second wing strip is equal to the horizontal tail semi -span. Given a modified variable-strip 
aerodynamic framework, model inaccuracy is introduced if consistent use of the second wing strip as a reference 
point when computing the downwash angle is employed, as the y-location of the outboard boundary of the second 
wing strip fluctuates with number of strips. Rather, a minimum difference approach is employed to identify the strip 
that retains an outer boundary y-value nearest to the semi-span of the horizontal tail. Computation of the downwash 
angle, a t i , is given as 

ds , 

(id 

where the sub-subscript r denotes the reference strip number previously discussed. 

In summary, while simple strip theory aerodynamics is applied to the gust response model, better utility is 
achieved by employing an aerodynamic model comprised of a user-defined number of strips. As a result of 
introducing a variable aerodynamic strip model, the force and moment equations are modified accordingly. 



Figure 7. Illustration of Enhanced Aerodynamic Strip Model 
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C. Extended Applications of Gust Inputs 

For completeness, a second continuous PSD function, given by Hugh Dryden, is added. Additionally, an 
extended application of the one-minus-cosine gust input is employed. The rational one-sided Dryden PSD function 
is expressed in the frequency-domain as 5 


* W (f) 
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(12) 


A comparison of the von Karman and Dryden PSDs shows that these spectra share similar shapes and have the same 
low-frequency asymptotes; however, at the high frequencies, the Dryden spectrum attenuates with a -2 slope, 
whereas the von Karman spectrum attenuates with a -5/3 slope. An extended application of the discrete gust input, a 
positive one-minus-cosine shape immediately followed by a negative one-minus-cosine shape, is employed in the 
present analysis and is computed in the frequency-domain as 
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The inverse Fourier transform of Eq. (13) yields a time-domain doublet one-minus-cosine gust input. The inclusion 
of an additional continuous frequency-domain input option and a doublet one -minus-cosine time-domain gust input 
enhances the gust response program by providing additional alternatives for analyzing aircraft responses to gusts. 


D. Modification of Horizontal Tail Mode Shapes and Capability to Animate Total Vehicle Deformation 
Response 

In Ref. 2, the single-strip horizontal tail is restricted to be rigid and is given pitch and plunge mode shapes equal 
to unity, respectively. The connection between the fuselage and the horizontal tail remains straight. 2 Given these 
restrictions, an animation of the total vehicle deformation response shows that the horizontal tail does not translate 
as a consequence of vehicle response in pitch and fuselage bending modes. In order to correct this response 
behavior, the pitch mode is defined to include the strips along the horizontal tail. Comparable to the modification of 
the horizontal tail pitch mode, a modified rear fuselage bending mode shape is applied to allow motion of the 
horizontal tail. In this modification, the horizontal tail is attached to the aft-most mass point along the fuselage. 

Analysis of the total vehicle deformation response to a gust input with a visual application is particularly useful. 
Developed in this analysis is the capability to produce such an application through gust response animations. 
Animation of the total vehicle deformation response to a gust input requires expansion of the existing structural 
model to include the z-direction displacements at the corner points of each strip of the aerodynamic model. 
Additionally, given the modification to include a user-defined number of strips, the expanded structural model 
requires the capability to compute the total deformation for a variable structural model. 

A stationary illustration of the capability to animate the total vehicle response to a gust is given in Fig. 8. In this 
depiction, the rigid body modes have been de-emphasized while the flexible modes have been exaggerated for 
visualization purposes. The horizontal tail has been omitted from this illustration. Frames are captured at t = 0.02 
seconds, t = 0.50 seconds, and t = 1 second. 


E. Addition of Aircraft Wing and Horizontal Tail Control Surface Capability 

While the application of active controls is beyond the scope of this research, the inclusion of a control surface 
capable of deflection is provided in the present analysis. Control surfaces are considered to be rigid and an infinitely 
stiff actuator is assumed. The implemented control surfaces are of user -defined span and chord and are positioned at 
user-defined wing and horizontal tail span locations. For computational simplification, this algorithm approximates 
the user-defined control surface position using a minimum difference approach so as to confine the inboard and 
outboard control surface boundaries to the boundary of the nearest strip. An illustration of a 20%-chord wing control 
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surface positioned with the outboard edge at 80% wing span and the inboard edge at 60% span for a ten-strip wing 
configuration is given in Fig. 9. 



Figure 8. Animation of Vehicle Deformation Response to Gust Input 



Figure 9. Wing Control Surface (10-Strip 
Configuration) 

IV. Results and Discussion 

A comparative evaluation is provided for analysis of the correlations and/or discrepancies between the revised 
gust load program and the original program. This comparative analysis serves as validation of proper 
implementation of these revisions. For example, a comparison of the time histories of the load responses obtained 
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from Ref. 2 and the time histories of the load responses obtained from the modified program, given a discrete one- 
minus-cosine gust input, is shown in Fig. 10. Fig. 10 shows that the baseline case is not disrupted by the revisions. 
Additional results illustrating the implementation of various enhancements can be found in Ref. 6. 

In conclusion, the following improvements and enhancements were successfully implemented, thereby 
furthering the cabilities presented in the original gust load analysis program: 1) aircraft wing and horizontal tail 
taper; 2) horizontal tail sweep; 3) variable number of strips; 4) dryden turbulence spectrum; 5) doublet one -minus- 
cosine gust input; 6) control surfaces; 7) animation of total vehicle deformation response; and 8) modified horizontal 
tail mode shapes. 



Figure 10. Comparison of Load Responses due to One-Minus- 


Cosine Gust Input 
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